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We have modeled transport properties of nanostructures using the Green's function method within 
the framework of the density-functional theory. The scheme is computationally demanding so that 
numerical methods have to be chosen carefully. A typical solution to the numerical burden is to 
use a special basis-function set, which is tailored to the problem in question, for example, the 
atomic orbital basis. In this paper we present our solution to the problem. We have used the finite 
element method (FEM) with a hierarchical high-order polynomial basis, the so-called p-elements. 
This method allows the discretation error to be controlled in a systematic way. The p-elements 
work so efficiently that they can be used to solve interesting nanosystems described by non-local 
pseudopotentials. 

We demonstrate the potential of the implementation with two different systems. As a test system 
a simple Na-atom chain between two leads is modeled and the results are compared with several 
previous calculations. Secondly, we consider a thin hafnium dioxide (Hf02) layer on a silicon surface 
as a model for a gate structure of the next generation of microelectronics. 

PACS numbers: 72.10.-d, 71.15.-m, 73.40.-c 



I. INTRODUCTION 

Using small nano-scale lithographic structures, atomic 
aggregates and even single molecules, it is possible to fab- 
ricate new kind of electronic devices [lj . The function and 
scale of these devices is based on quantum-mechanical 
phenomena, and cannot be described within the classical 
regime. Of particular relevance are the electron trans- 
port properties of these nanoscale devices, as this will 
determine their effectiveness in, for example, a new gen- 
eration of transistors. As the experimental work on these 
devices grows, increasing emphasis is placed on develop- 
ing a matching theoretical description 0, Q . Although 
some efforts have included a full description of an elec- 
tronic circuit 0, H , current research is mainly focused to 
study single electronic components. 

Density-functional theory (DFT) is widely used in 
atomistic modeling of materials properties and recently 
also properties and phenomena in nanostructures. The 
power of DFT is in its capacity to treat accurately sys- 
tems with a hundreds of atoms, yet retain a full quantum- 
mechanical treatment. Although the full justification of 
use of the DFT in electron transport calculations is de- 
bated 0, Q we adopt it as a practical scheme to describe 
the real systems and devices. 

In the Kohn-Sham scheme of DFT the electron den- 
sity is calculated using single-particle wave functions. 
The explicit use of the wave functions in constructing 
the density suffices well in two kinds of systems. Either 
the system has a repeating structure so that it can be 
modeled with periodic boundary conditions or the sys- 
tem is so small that it can be calculated as a whole. In 
nano electronics, however, a system consists usually of a 
small finite part, the nanostructure, which is connected 
to the surrounding infinite leads. If one enforces periodic 
boundary conditions even a large repeating super cell or 



calculation volume can cause finite-size effects with spu- 
rious results for electron transport. 

A commonly used solution to this problem, which we 
have also employed, is to combine DFT with the Green's 
function formalism The Green's functions are first 
constructed for the semi-infinite leads by using the ana- 
lytically known or easily calculated wave functions. Once 
the Green's function for the combined nanostructure and 
leads is constructed, the wave functions are no longer 
needed explicitly. This makes it possible to use open 
boundary conditions between the nanostructure and the 
lead. In this way we have an effectively infinite system 
without periodicity, making the finite-size effects small. 
It is also possible to calculate the electric current through 
the system for a finite bias voltage between the leads in a 
self-consistent manner with the electron density. The en- 
suing model for the current is analogous to the Landauer- 
Biitteker model @. We have used non-local pseudopo- 
tentials for modeling atoms, and the ideal metal "jellium" 
model for the leads. The charge density in the leads can 
be varied according to the conducting properties of the 
leads we wish to model. 

The use of Green's functions instead of the explicit use 
of wave functions is computationally demanding. This 
is why special care has to be taken in choosing the nu- 
merical methods. The first implementations used tight- 
binding methods 0, but a more typical solution is 
to expand the Green's functions in a special basis tai- 
lored for the system. Common examples are localized 
atomic orbitals El E|, an 0(N) optimized basis [3, 
a wavelet basis |13|, full-potential linearized augmented 
plane- waves |l4j . maximally localized Wannier functions 
|15|. a finite-difference method ^(|, and linear a finite- 
element method EJ- O ur solution is to use the finite- 
element method (FEM). It allows a systematic error con- 
trol which is especially important in transport problems 
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as there are many different properties which must be 
monitored. For example, the pole of the Green's func- 
tion can cause numerical problems. According to our 
experience electronic tunneling in particular is sensitive 
to numerical accuracy. 

Besides systematic error control, the FEM has also 
other good properties which makes it a natural method 
for transport problems. It is a flexible method which 
allows one to take into account the geometry of the 
nano device exactly. Special boundary conditions are 
easy to derive without mixing the model with the nu- 
merical method and their implementation is straightfor- 
ward. Moreover, the local nature of the basis produces 
sparse matrices for which efficient solving methods exist. 
Varying the size of the elements can be used to reduce 
the number of the basis functions and, consequently, the 
size of the system as compared to the finite-difference 
method. This is especially true for the high-order p- 
method. Finally, there exists a lot of theoretical work 
together with tested and reliable tools, such as mesh gen- 
erators and optimized linear solvers. These are used as 
standard building blocks for any FEM implementation 
granting easy access to state-of-the-art algorithms. Us- 
ing the FEM new theoretical or numerical ideas are easy 
to implement and test. 

The structure of the paper is as follows: in Sec. 1TT1 we 
describe the model itself in detail, including a discussion 
of the formalism of our implementation; in Sec. IIIII we 
apply the model to two example systems, a Na atom 
chain and Hf02-Si interface between two leads. In Sec. 
IIVI we summarize the work. In this paper we use atomic 
units in all equations. 



II. MODEL 

The schematic picture of our model is shown in Fig.^ 
Actually, the figure present our second test case, the 
Hf02-Si interface between two leads. We have an atom- 
istic nanostructure between two semi-infinite leads. The 
system is divided into three parts, fl being the calcula- 
tion volume, and il^ and flu are left and right leads, 
respectively. The boundaries d£l L / R are open so that 
electrons can penetrate through them without any re- 
flection or refraction. We use the DFT to model electron 
interactions. The basic quantity, the electron density, is 
calculated from single-particle Green's functions. Then 
we use the density to calculate the effective potential as 

V eS = V ext + V C + V xc + Vbias + Vnh (1) 

where Voxt is the external potential caused by positive 
background charges, local parts of the pseudopotential 
operators and the potential outside potential barriers. 
V c is the Coulomb Hartree interaction part, and V xc is 
the exchange-correlation part which we calculate using 
the local-density approximation parametrized by Perdew 
and Zunger [jjl ■ Vt,i as sets the boundary conditions 




FIG. 1: Schematic picture of the model. The Hf02 interface 
is used as an example. The small and large gray gray spheres 
denote the Hf and Si atoms, respectively, and dark spheres the 
O atoms. The gray volumes are the jellium leads. The system 
consist the volumes Ql, f2, and Hp and of the boundaries 
OQl, dflp, dtlpi, d£lp2, dflp3, and dQp4- 



if a bias voltage is applied. V n \ is the nonlocal part of the 
pseudopotential operators. 

The Hartree potential is calculated from the modified 
Poisson equation 

V 2 V: - k&V* = 4tt( P+ - p) - k 2 P Vr\ (2) 

where fcp is an adjustable parameter, kp does not af- 
fect the final self-consistent result, but the stability and 
convergence of iterations are improved pfl | , because the 
Coulomb potential due to charge redistribution between 
adjacent iterations is screened. The non-local pseudopo- 
tential is an operator given by 

V nl v(r) = yV m Cz, m (r) / Q, m (r')v(v')dr', (3) 
l, m Jn 

where e; m and Cz,m( r ) are defined using the Troullier- 
Martins pseudopotentials IMSl- Eq. © uses the pro- 
jection of the function v(r) (arbitrary function, which in 
practical calculations is a basis function) on the atomic- 
specific function C,i >m depending on the quantum numbers 
I and m corresponding to the angular momentum. 

We have implemented the guaranteed-reduction-Pulay 
[23| method for the mixing of the self-consistent itera- 
tions. It uses potentials from the five previous iterations 
for computing a new potential in such a way that the pre- 
dicted norm of the potential residue is minimized. The 
simplest mixing scheme in which potentials are mixed 
with a linear feed-back coefficient does not work well in 
open systems. The calculations are rather unstable so 
that quite a small feed-back coefficient has to be used. 
This is because the net charge in the calculation volume 
f2 varies during the calculations. 
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A. Green's function model 

The details of the Green's function method for electron 
transport in nanostructures are explained, for example, 
in Ref . 24] . Here we give only a short introduction to the 
equations to be solved. The retarded Green's function G r 
is solved from the equation 



(c-ff(r))G' r (r,r';c)=5(r-r'), 



(4) 



where u> is the electron energy and H is the DFT Hamil- 
tonian of the system, 



H(r) = -\V 2 + V ctf (r). 



(5) 



When we know G r we can calculate the so-called lesser 
Green's function G < . In the equilibrium when no bias 
voltage is applied over the nanostructure it is obtained 
from 

G<(r, r'; w) = 2f L/R (u) G r (r, r'; cj), (6) 

where /l/r are the Fermi functions of the leads. In the 
equilibrium, Jl = f R . For a finite bias voltage Jl/r 
we take into account the bias and a more complicated 
equation for G < has to be used. To obtain it we write 
Eq. J2J| in the form 

(u>-H - ££(u;) - £fc(w)) G r (r, r'; u) = S(r - r'), (7) 

where Hq is the Hamiltonian of the isolated volume £1 
and ^ r L / R are the so-called self-energies of the leads. We 
also define 

iV L/R = ^ L/R - EJ, R = 2z Im(E£ /fl ), (8) 
and can solve G < for a finite bias voltage as 
G<(r,r';w) = 

-if R (u) I I G r (r 7 VR;iu)TR(rR,r' R ;u;) 

xG a (r R ,r';u)dr R dr' R (9) 
-if L (u>) [ I G r (r,r L ;u)T L (r L ,r' L ;cj) 

JdU L JdCl L 

x G a (r' L ,r';uj)dr L dr' L . 

The first and second terms correspond to electrons orig- 
inating from the right- and left leads, respectively. The 
electron density is calculated from 

-1 r°° 

p(r) = — Im(G<(r,r; W ))da; (10) 

and the tunneling probability from 

T(u>)=[ Iff r L (r L ,r' L ;uj)G r (r' L ,T R ;uj) 
Jda L Jon L Jon R Jan R 



xT R {r R , r' R ; w) G a (r' R , r L ;uj) dr L dr' L dr R dr' R , 



from the values of the functions at the boundaries dfl^ i R . 
Finally the current is determined as 



I=- I T(u)(f L (uj)-f R (uj))duj. (12) 

7T ' 



We use the FEM in the numerical implementation. 
Therefore we first cast Eq. (@J in the variational form 
with open boundary conditions (for the derivation, see 
Ref. 2S]) as 



- Vv(r) • -VG r (r,r';w) 

+ v(r) [u - V eS (r)] G r (r, r'; w)} dr ( 13 ) 

- < ± L G r ,v > - < ± R G r ,v> 
■ v(r'), 



where the self energy-operators 



< ErG r ,t; >= 



d 2 ge(r' L ,r L ;w) 
dxiLdn'r 



(14) 



v(r L ) dv' L dv L . 



(11) 



Above, g e is the Green's function of the semi-infinite lead 
in the domain / R with the zero- value condition on the 
boundary dflL/R. In our implementation the leads are 
described by a uniform positive background charge and 
therefore g e can be calculated partly analytically. Thus 
our model means that the leads are of some kind of ideal 
generic metals. The important interface between the 
nanostructure, e.g. a molecule, and the actual metal- 
lic lead can be described accurately by including some 
lead metal atoms in the computational domain Q. It is 
also possible to use fully atomistic leads by calculating 
numerically g e for them. 

Note that the Eqs. I|13|) and i|14|) are analogous to those 
derivations of the open boundary conditions in which 
truncated matrices [l| are used. In the continuum limit 
these two forms give the same results. However, the weak 
form is more natural in the FEM formulation and more 
suitable for theoretical purposes when analyzing nonlin- 
ear partial differential equations. It is also straightfor- 
ward to use and the error control is systematic. Note 
that this formulation can be used with any continuous 
basis set, not only with the FEM. In the context of basis 
set methods, the weak formulation case is known as the 
Galerkin method. In practice the Green's functions are 
approximated with respect to this basis so that 

N 

G r (r,r';cj) a £ ff«(w)&(r)&(r'). (15) 

»>i=i 

The coefficients gij(uj) can be solved from (fTSfl by choos- 
ing v = 4>k and evaluating the equations. 
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B. Finite-element p-basis 

In the FEM we partition the calculation volume to (in 
our case, tetrahedral) sub-domains called elements and 
the basis functions fa are constructed using globally con- 
tinuous (but not necessarily continuously differentiable) 
piecewise polynomials with respect to the finite element 
mesh. This gives both unique flexibility of the approxi- 
mating functions as well as completeness of the basis with 
respect to almost any norm. Each basis function has a 
support that is concentrated to only a few neighboring 
elements. This makes the basis local and results in sparse 
system matrices. 

There are several options of how to choose the finite- 
element basis and one has to be careful in achieving ac- 
ceptable accuracy. The simplest basis is the linear one. It 
is easy to implement and works well, especially in systems 
with rapidly varying functions. A typical improvement 
to this basis is to use node-based higher-order elements. 
These elements converge faster to a smooth solution than 
the linear ones. However, practically only relatively low 
orders, two and three, can be used because of numerical 
stability problems. 

In this work we have used so-called hierarchical p- 
elements. They also span higher-order polynomials, but 
the choice of the local basis ensures that stability prob- 
lems do not appear. This is because the basis functions 
are chosen so that their derivatives are close to orthogo- 
nal in the L2-norm. The hierarchical nature also makes 
it easy to change the order of the basis from element to 
element within the same mesh. 

The actual FEM implementation consists of a refer- 
ence element and reference basis that are mapped sep- 
arately to each of the elements of the mesh. Our ref- 
erence element is a tetrahedron with nodes at the co- 
ordinates 1 : (-1,0,0), 2 : (1,0,0), 3 : (0,\/3,0) and 

4 : (0, -^,2^/|). One can easily show that there exists 
an affine map taking the reference element to any of the 
tetrahedron. The order of our basis is p meaning that in 
each element polynomials of the order of p are employed. 
The basis is constructed hierarchically. First, there are 
four linear node basis functions inside the elements shar- 
ing a common node. In the reference element they are 



L ' = - [ "-js 




functions. E.g. for the edge between the nodes 1 and 2 
A^ = LiL 2 (pi(L 2 - Li), i = 2, ...p. (17) 
Here one usually sets 

(18) 



i(0 



2(2i - 1 



Above Pi is the Legendre polynomial of the order of i. 
Third, we have 2(p— l)(p— 2) face functions. For example 
for the face between the nodes 1,2, and 3 they are 



N, 



(1,2,3) 



L 1 £ 2 i 3 i , <(i2-£i)Pi(2L 3 -l) > 



where £ fi and £ are the cartesian coordinates of the ref- 
erence element. Secondly, for p >1 we have 6(p — 1) edge 



tj j ~±— f H-i o ~n (19) 

i,3 = 0, ...,p-3, i+j = 0, ...,p-3. 

Fourth, we have s(p — l)(p — 2) (p — 3) bubble functions, 
which are supported only in a single element each. These 
are 

N iJtk = L 1 L 2 L 3 L 4 P l (L 2 - Li)x 

P j (2L 3 -l)P k (2L i -l) (20) 
k = 0, ...,p-4, i + j + k = 0, ...,p- 4. 

When p-elements are used one must take care of the 

continuity of the basis. This is because, for example, 

(l 2) 

the local basis function AQ ' has an orientation on the 
boundary. The basis includes the function tp(L 2 — L±), 
not (p(L± — L 2 ), which would be another possibility. This 
means that all the edges in the mesh have to have in- 
formation about the direction. Otherwise there is very 
likely a continuity problem on some boundaries. In prac- 
tice, for tetrahedral elements the orientation problem can 
be handled for arbitrary finite-element meshes using only 
two reference elements |26| . 

The benefits of selecting the basis described above are 
rather clear. The polynomial basis is very easy to real- 
ize and has good approximating properties. For smooth 
solutions the p-basis is known to give exponential conver- 
gence rates with respect to the number of basis functions 
used. In the DFT methods the theory is typically devel- 
oped to the direction that the solutions are as smooth 
as possible. For example, pseudopotential operators are 
designed so that they produce as smooth an electron po- 
tential as smooth as possible. This is because the plane- 
wave basis set needs smooth solutions in order to work 
efficiently. On the other hand, in the case of non-smooth 
solutions one can benefit from the piecewise nature of the 
FEM basis allowing one to approximate even singular so- 
lutions to some extent. Moreover, the finite element mesh 
can be refined in regions where solution changes rapidly. 
When modeling molecules there is also a lot of empty 
space in the calculation domain. It is then practical to 
use large elements in the empty space and smaller ones 
near atoms. 
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C. Linear algebra methods 

The use of the Green's function method is compu- 
tationally demanding in comparison to explicit wave- 
function methods. Since the main computational burden 
of our method is to find a subset of the coefficients of 
the Green's function in question, a special consideration 
must be taken when choosing the methods of linear al- 
gebra to be used. The eigenvalue problems in explicit 
wave-function methods are typically solved by iterative 
methods. In our case it is better to use direct solvers, 
because a set of linear equations needs to be solved. We 
have opted for the frontal method widely used in the 
solution of sparse linear systems 27, 28] and extremely 
suitable for finite-element matrices. The actual imple- 
mentation is ME47 of the Harwell Subroutine Library 
(HSL) [H (see Refs. El H3 for other similar ap- 
proaches). In the frontal method, one first finds a per- 
mutation of the sparse matrix aiming to minimize the 
fill-in resulting from the factorization process. Next, a 
LU-decomposition (or Cholesky-decomposition for sym- 
metric problems) of the matrix A is found, and finally, 
two systems with triangular coefficient matrices, Lz = b 
and Ux — z (where U ~ L T for symmetric problems) are 
solved. To find all the required coefficients of the solution 
we must vary the right-hand side b of the system. 

For a three-dimensional problem the size of the linear 
system can grow so large that the CPU-time and memory 
requirements of different systems have to be addressed. 
The main question is how large systems can be calculated 
using these methods so that the calculation time for a 
single self-consistent iteration is not too large. Currently 
a system of several tens of thousands of unknowns can 
be solved in a commodity-CPU cluster environment. 

In detail, the Green's function method includes a com- 
putation of the elements for the inverse of a sparse ma- 
trix, so that the calculation time requirements increase 
relatively fast with the system size. A classical complex- 
ity result for the solution (and inversion) of a general 
N x N system with a direct method is 0(N 3 ). How- 
ever, for sparse systems and modern frontal methods this 
bound is too pessimistic |3^|. The CPU time require- 
ment depends on the fill-in of the inversion problem. For 
very simple cases one can show that the key statistic 
of the problem, the number of non-zeros (nnz) present 
in the factors L and U, satisfies nnz(L) ~ nnz(U) = 
0(N log (AT)) |23. Then the solution of each of the 
systems requires 0(nnz(L) + nnz(R)) floating-point op- 
erations, and in the worst case we must solve these 
with N different right-hand sides effectively giving us 
the inverse of the matrix A, so that the total cost is 
0(N(nnz(L) + nnz(R))) . However, in modern computer 
systems the complexity is not the only relevant measure 
since the performance may be highly nonlinear (see, e.g. 
for an example on BLAS-tuning). 

Another topic related directly to the performance of 
modern computer systems is the relation between proces- 
sor power and memory bandwidth. This is especially true 



for the computation of the Green's function where the ac- 
tual bottleneck is the lack of available memory bandwidth 
in commodity-based cluster systems used in calculations, 
not the floating-point performance of the processor itself. 

It is likely that a better performance can be achieved 
by upgrading several parts of the algorithms. First, the 
current parallel solver is implemented using the Mes- 
sage Passing Interface (MPI) j3g. However, in Symmet- 
ric multiprocessor (SMP) systems it is likely that well- 
designed OpenMP [33 (or similar) parallelism would re- 
duce the need for data transfer and thus increase perfor- 
mance. It would also decrease the memory requirements 
of the problems. Second, at the moment the solution 
of the Green's function is computed varying one vector 
on the right-hand side at a time. A better performance 
could be obtained if the equations could be solved for 
multiple right-hand sides at a time allowing the use of 
BLAS3 routines. Finally, it is likely that computations 
would benefit from a computer system having a larger 
memory bandwidth than our present commodity-based 
one. 



III. EXAMPLE SYSTEMS 

A. Atomic wire 

Using the atomic force microscope or the mechanically 
controlled break junction technique, a chain of atoms can 
be made of certain metals j^. It has been observed 
that the conductances of atom chains vary as a function 
of the number of atoms in the chain 39]. The conduc- 
tances of these systems have been studied also theoreti- 
cally in several works. In order to benchmark our results 
against other calculations, we use Na-atom chains as test 
systems. They have been simulated in several previous 
studies |4(J, HH H3, |43j using different models. According 
to these calculations the conductances of the wires show 
even-odd oscillation as a function of the number of atoms 
in the wire. 

In our setup, the Na-atom chain is located between two 
leads, with the lead shape defined by a 70° cone angle 
(see Fig.[2J|. We consider two different connections of the 
atom chain to the electrodes. In model A we have just 
three Na atoms between the jellium leads. This resembles 
closely the system used in Ref. ^l|. In model B, there are 
four Na atoms at the tips of the leads in a square form. 
This makes the connection between the atom chain and 
the leads more realistic. This kind of structure is modeled 
also in Ref. |42j. 

The conductances as a function of the number of chain 
atoms for systems A and B are shown in Fig. In the 
Na-atom chain, electrons have only one conducting mode 
so that the conductance can be one conductance quan- 
tum 2e 2 /h at maximum. Both systems A and B exhibit 
conductance oscillations as a function of the number of 
atoms. These oscillations arise from resonance states in 
the atom chain. Depending on the position of the reso- 
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FIG. 2: Two models of the Na-atom chain, a) In model A, 
Na atoms are directly connected to the cone-shape leads, b) 
In model B, there are four Na atoms as squares at the tips of 
the leads. 
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FIG. 3: Conductance through the Na-atom chain as func- 
tion of the number of Na atoms in the chain. The results 
for systems A and B (see Fig. |5J are denoted by circles and 
stars, respectively. For system A with three Na atoms, results 
corresponding to 50° and 90° cone angles a are also given. 



nances relative to Fermi-level the conductance has either 
a maximum or minimum value, so that the maxima and 
minima correspond to approximately half and fully oc- 
cupied resonance states, respectively. The oscillation is 
within the range of 0.9 - 1.0 x2e 2 /h for system A and 
0.6 - 1.0 x2e 2 /h for system B. The difference between 
the oscillation amplitudes is due to different strengths of 
the connection of the chains to the leads. System B has 
a weaker coupling to the leads than system A. Weak con- 
nections make the resonances also sharper, as is seen in 
the tunneling probability in Fig. In contrast to Ref. 
[44l | , we do not see a strong lead-shape dependence in the 
conductance. The widening of the cone angle lowers the 
conductance as the edges of the wire become sharper. 
The electron tunneling probabilities through chains 
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FIG. 4: Tunneling probability from Eq. (I I 1 1 for three- and 
four- atom chains between two semi-infinite jellium leads. The 
solid and dashed lines correspond to systems A and B shown 
in Fig. |5J respectively. The Fermi-level is marked by dashed 
vertical lines. The cone angle a=70° 

of three- and four- atom systems A and B are shown 
in Fig. 21 The probability function T{u) is defined in 
Eq. i|ll|) ■ The conductance of the system in the zero-bias 
limit can be read at the Fermi-level. Here, as well as in 
Fig |3 we see that the conductance oscillations for sys- 
tems A and B are in a different phase. This is because 
in system B the atom chain is effectively shorter than in 
A, as the first and the last chain atom are partly inside 
the square of the four Na atoms. 

When we compare the conductance oscillations of sys- 
tem A (see Fig. |2J to those in Ref. ^l| obtained by us- 
ing semi-infinite jellium leads with planar surfaces (a — 
180°), we see that the even-odd oscillations in the con- 
ductance are in the same phase. In the case of system 
B we can directly compare the tunneling probability of 
Fig. H with those in Ref. H3 where the atom chain is con- 
nected also through a square of four Na atoms to jellium. 
The phase and the amplitude of the conductance oscil- 
lations of these results are in good agreement with our 
values in Fig.|3J Now that we have satisfied ourselves that 
the method provides a good model for electron transport 
we can consider a more interesting and demanding ex- 
ample. 



B. Thin insulating layer 

The general increase in the performance of microelec- 
tronic devices in the past few decades has been made 
possible by continuous transistor scaling - based on a re- 
duction in the thickness of the gate dielectric in typical 
metal-oxide-semiconductor field-effect transistors (MOS- 
FET). At present the process has reached a bottleneck, 
as further reduction leads to a large increase in leakage 



current due to direct tunneling across the thin silicon 
dioxide (Si0 2 ) layer. Several pos sible approaches to re- 
solve this are being considered [43 , but retaining conven- 
tional MOSFET design remains an economically attrac- 
tive choice, and a leading option is just to replace SiC>2 
with another oxide of higher dielectric constant (high-fc). 
A high-fc oxide would provide higher effective capacitance 
to a comparable Si02 layer, hence allowing thicker layers 
to be used to reduce losses due to tunneling. The specific 
choice of oxide is determined by a set of requirements |4rj| | 
based on both the intrinsic properties of the grown oxide 
and its integration into the fabrication process, and at 
present hafnium oxide (Hf02) remains a leading candi- 
date. 

In order to study to transport properties of thin Hf02 
films we have simulated the growth of the oxide on a sil- 
icon surface via first principles molecular dynamics |47j . 
Here we consider three model interfaces: (i) a nonstoi- 
chiometric oxide interface (C), which is basically metallic 
due to Hf-Hf and Hf-Si bonds across the interface; (ii) a 
stoichiometric oxide interface ( D), which has a localized 
state in the band gap due to a few Hf-Hf bonds; (iii) a 
more idealistic interface (E), which remains insulating if 
no defects are present. The last model is based on the 
interface used in Ref . 0] , but slightly reduced in size to 
make it computationally manageable 49j . These models 
were calculated with periodic boundary conditions with 
k-points on the boundaries cW^pi/p2/P3/P4- The effective 
potentials have been calculated for systems C and D us- 
ing the gamma point, and for system E, four fc-points. 
All the tunneling probabilities T(ui) are calculated using 
four fc-points, which was enough to converge the proba- 
bilities to a good accuracy. 

As shown in Fig. ^ the interface models are positioned 
between two leads. The charge density in the leads is 
chosen so that in the right lead r s = 2 (electron density 
n e — 3/(47rrf)), representing a metal, and in the left one 
r s = 3.1, representing doped-silicon - as in a standard 
MOSFET design. 

The tunneling profiles of the systems are shown in 
Fig- El Here it is seen that systems C and D show clearly 
metallic behavior, with a large tunneling probability at 
the Fermi energy. Although in principle, the stoichio- 
metric interface (D) has a much lower density of metallic 
bonds, it is clear that both in interfaces C and D around 
two channels dominate the transport. The localized de- 
fect state in the band gap of system D plays an equivalent 
role in transport to the metallic bonds in interface C. 

As expected, the tunneling probability for the more 
ideal interface E is an order of magnitude smaller at the 
Fermi energy than those for interfaces C and D. Yet we 
also see that it remains significant - this is largely due 
to the structure of the interface itself |4g. Although 
bulk HfC>2 is a wide-bandgap insulator, at the interface 
it exists as almost tetragonal HfSiC>4, and the effective 
band gap is actually smaller than that of bulk silicon 
below the interface. This means that there is a negative 
conduction band offset between silicon and HfC>2 , and no 
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FIG. 5: Tunneling probability T(lj) through thin Hf02 lay- 
ers. Results for system C (dotted line) - nonstoichiometric 
interface, system D (solid line) - stoichiometric interface and 
system E - ideal interface are shown. The inset shows T(ui) 
normalized with the conducting area enabling the comparison 
of actual insulating properties of different systems. 

real barrier for leakage. Although some of this is caused 
by the underestimation of the band gap in the DFT, this 
also reduces the silicon band gap (although the effect is 
not systematic). 

0.3i . . . . . 1 




Position along transport direction (a o ) 

FIG. 6: Change of the average effective potential in inter- 
face D when a 0.25 V bias voltage is applied over the Hf02 
layer. In every position along the transport direction the effec- 
tive potential is averaged over the perpendicular-coordinates. 
Atom positions are indicated: open circles are silicon, filled 
circles, and stars hafnium and oxygen atoms, respectively. 
The gray areas mark the positions of the leads. 

The poor performance of interface D can also be seen 
in its capacity for dropping the potential. Fig. shows 
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the potential change for 0.25 V applied bias voltage. The 
potential drop across Hf0 2 is less than 0.05 eV, demon- 
strating that the oxide hardly perturbs the electron flow 
from the right lead. The potential drops fastest at the 
right hand side of HfC>2-layer where pure HTO2 exists, 
and much more slowly in the thin layer of SiC>2 formed 
due to diffusion of oxygen. The large drop at the lead 
and silicon atoms is just an artifact of the boundary con- 
ditions of the Coulomb part of the effective potential. 

In the rigid band approximation (used for example in 
Ref. H3) it is assumed that the shape of the tunneling 
probability stays constant and is only shifted in energy 
so that T(ui, Vbias) = T(u> + i]Vbm S ), where 77 is the ratio 
of potential drop at the other end of the nanostructure to 
the total drop over the nanostructure. In Fig. we have 
studied how well this approximation works for interface 
D. The curves are plotted so that the zero-bias Fermi 
level is in the middle of the left and right Fermi levels of 
the biased interface. This corresponds to the symmetric 
case with 77 = 0.5. We see that the tunneling proba- 
bility curves roughly coincide. This indicates that po- 
tential drops symmetrically over the nanostructure and 
the rigid-band approximation gives a rather reasonable 
result. 
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FIG. 7: Tunneling probability through a thin Hf02 layer. 
Results for interface D with 0.25 V bias (solid line) and zero 
bias (broken line) voltages shown. The vertical lines show the 
positions of the Fermi levels in the left and right leads for 
the biased system. Between them is located the "so-called" 
conductance window [see equation l|12|l ]. The Fermi level of 
the non-biased system is in the middle of these lines. 

The above results show, in agreement with previous 



calculations [48( that tunneling through a more ideal, 
insulating interface is still significant due to a negative 
band offset with silicon. Since the only Hf02 interfaces 
providing significant band offsets to silicon were built 
very idcalistically (i.e. assuming no significant atom mi- 
gration nor interfacial Si02 growth) |5fl l5l||. this indi- 
cates that fabricating a good interface directly between 
silicon and Hf02 is very difficult. A more viable alterna- 
tive maybe to sacrifice somewhat in dielectric constant, 
and grow Hf02 onto a pre-existing SiC>2 layer. These 
possibilities will be explored in more detail in a further 
work |47j . 

IV. CONCLUSIONS 

In this paper we present a finite-element implemen- 
tation of the non-equilibrium Green's function method 
which is combined to the density-functional theory. Al- 
though the Green's function method is computationally 
demanding, we demonstrate that by using hierarchical 
p-elements, large, physically relevant systems become 
tractable. More importantly, our method offers a much 
more rigorous control of accuracy than is usually possible 
in transport calculations. 

We demonstrate the functionality of our implementa- 
tion with two kinds of systems, the sodium atom chain 
wire and the silicon-Hf02 interface. For the atom chain, 
we show that the method reproduces the previous re- 
sults of other Green's function transport methods. This 
gives us confidence to apply it to the more complex sys- 
tem: a thin layer of hafnium oxide on a silicon substrate. 
Here we show that the transport properties are an even 
more sensitive indicator of the role of defects than the 
electronic structure. Comparison of stoichiometric and 
non-stoichiometric HTO2 oxide layers demonstrates that 
even one or two defects in a stoichiometric interface can 
result in tunneling comparable to that of a fully metallic 
non-stoichiometric interface. 
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